---
title: "R13.Expropriation_IV"
---

# Load Packages
```{r}
library(tidyverse)
library(estimatr)
```

# Load Data
```{r}
load("Data/data_final.Rdata")
load("Data/gnr_panel.Rdata")
load("Data/ideo_list.Rdata")
```

# Set Covariates
```{r}
covs_list <- c("pr_ha_ord", "as.factor(DistN)", "ln_PopDens_70", "Cities_10km", "ln_Confs_All", "PrRecH_Adm3", "Isolados1950", "PrWheatAr_1972", "LandIneqRat", "TerrainSlopeIndex", "MedAltitude", "Communist_1975", "MinDistPCP_Map", "log_gnr_cumul_23_Z")
```

```{r}
elections <- c("AC75", "AR76", "AR79", "AR80", "AR83", "AR85", "AR87")
```

# Expropriation as IV
```{r}
  
baseline_list <- list()  
for (i in 2:7){  # for elections between 1976 and 1987

aggregated_1 <- ideo_list[[1]] %>% select(-election)  # Get vote shares for baseline year (1975)
colnames(aggregated_1)[3] <- paste(colnames(aggregated_1)[3], "1", sep = "_")
 
aggregated_2 <- ideo_list[[i]] %>% select(-election)  # Get vote shares for year of interest
colnames(aggregated_2)[3] <- paste(colnames(aggregated_2)[3], "2", sep = "_")

# First Difference
aggregated_fd <- aggregated_2 %>% inner_join(aggregated_1) %>% rowwise() %>%# Subtract 1975 vote shares from vote shares of year of interest
  mutate(FD = votes_2 - votes_1) %>%
  select(DICOFRE_Z, ideology2, FD) %>%
  pivot_wider(c(DICOFRE_Z), names_from = "ideology2", values_from = "FD")

data_sample_g <- data_final %>% left_join(aggregated_fd) %>% # Join vote share first difference and cumulative GNR events up to election year to main data set
  left_join(gnr_panel %>% filter(elec_number == i)) 

output <- list()  # Subset for (3) groups (land reform intensity)
for (g in 0:2){ 
 
output[[g+1]] <- data.frame(election = elections[i], group = g,  model = colnames(aggregated_fd)[-c(1)], pr_ha_ord = NA, se = NA, n = NA) %>% filter(model != "NA")

      # Run models for each ideology category (4)
      for (r in 1:nrow(output[[g+1]])){ # r<-1
      m.1 <- lm_robust(as.formula(paste(output[[g+1]]$model[r], paste(c(covs_list), collapse= " + "), sep = " ~ ")),
                data = data_sample_g %>% filter(exp_cat2 == g), weight = weights, cluster = subclass)
      
           output[[g+1]][r, "pr_ha_ord"] <- m.1$coefficients["pr_ha_ord"]
           output[[g+1]][r, "se"] <- m.1$std.error["pr_ha_ord"]
           output[[g+1]][r, "n"] <- length(m.1$fitted.values)
      }

}

output_list <- do.call(rbind, output)
baseline_list[[i]] <- output_list

}

# Plot Results

baseline_list_all <- do.call(rbind, baseline_list) %>% 
  mutate(model = ifelse(model == "Center", "Center-Left", model), 
         model = ifelse(model == "Left", "Far Left", model)) %>% 
  filter(group != 0)

years <- data.frame(election = c("AR76", "AR79", "AR80", "AR83", "AR85", "AR87"), 
           year = as.factor(c(1976, 1979, 1980, 1983, 1985, 1987)))
baseline_list_all <- baseline_list_all %>% left_join(years)
baseline_list_all$model <- factor(baseline_list_all$model, level = c("Communist", "Far Left", "Center-Left", "Right"))


ggplot(baseline_list_all, aes(x = year, col = as.factor(group))) +
  geom_point(aes(y = pr_ha_ord, shape = as.factor(group)), position = position_dodge(width = 1)) +
  geom_linerange(aes(ymin=pr_ha_ord - (qt(0.975, n)*se), ymax=pr_ha_ord + (qt(0.975, n)*se)), position = position_dodge(width = 1), size = 0.5) +
  geom_linerange(aes(ymin=pr_ha_ord - (qt(0.95, n)*se), ymax=pr_ha_ord + (qt(0.95, n)*se)), position = position_dodge(width = 1), size = 0.75) +
  geom_hline(aes(yintercept = 0), col = "red") +
  facet_wrap(~model, scales = "free_y") +
  xlab("Election Year") + ylab("Effect of Expropriation") +
  scale_color_manual(name = "Land Reform Intensity", labels = c("LR < 50%, or adjacent", "LR >= 50%"),
                     values = c("#636363", "#969696"))+
  scale_shape_manual(name = "Land Reform Intensity", labels = c("LR < 50%, or adjacent", "LR >= 50%"),
                       values = c(17, 15)) +
  theme_light() +
  theme(axis.text.x = element_text(angle = 90))

ggsave(plot = last_plot(),
             "Figures/Afg4d_ExpropReform.pdf",
             width = 9, height = 5)
  
```


